The goal of seismic processing is to image subsurface reflectors. During a seismic survey, a large amount of the seismic energy generated by a source does not travel vertically as body waves in the earth to reflectors and then return to the geophones, but rather travels horizontally through the shallow near surface of the earth from the source to the geophones. These strong surface waves, sometimes called ground roll, can mask the weaker reflections and should be removed prior to imaging. In other words, the waves of interest are the deeper body waves from reflectors, such as hydrocarbon reservoirs, and the body waves are obscured by overlapping surface waves.
Prior surface-wave removal methods attempt to exploit differences between body waves and surface waves. Two of the most apparent differences are in their relative velocities and frequencies. When multi-component data is available, differences in phase between two or more co-located seismic receivers can also be exploited to mitigate surface waves.
A problem with surface-wave mitigation methods is that some reflector energy is often removed along with the surface-wave energy. The removal may occur for several reasons, including incomplete separation between the velocities and frequencies of surface waves and body waves. Another reason is aliasing of the surface waves, which makes it difficult to remove them by classical Fourier (or dip separation) methods.
In Fourier or dip-separation methods, data are transformed into the Fourier (f-k or frequency-wavenumber) or Radon (τ-p or time-slowness) domain. In these domains, data having different dips in seismic records appear at different locations in the transformed domain. Fourier methods are classical and widely used, for example see O. Yilmaz, Seismic Data Processing, Society of Exploration Geophysicists, 1987, which is hereby incorporated herein by reference. Improvements on simple f-k methods include applying high-resolution spectral techniques to improve the distinction between waves of similar dips, see Iranpour, U.S. Pat. No. 6,834,236 B2, which is hereby incorporated herein by reference. Most such methods stem from the original work of J. Capon, 1969, “High-resolution frequency-wavenumber spectrum analysis,” Proceedings of the IEEE, vol. 57, no. 8, 1408-1418, which is hereby incorporated herein by reference.
Radon methods are a more recently developed transform and more computationally expensive, but are sometimes used for surface-wave filtering, for example see Robinson, U.S. Pat. No. 7,239,578 B2, which is hereby incorporated herein by reference. One advantage of the Radon transform is that it does not require regularly spaced trace sampling in the seismic record. Fourier methods suffer from irregular sampling in addition to the previously mentioned aliasing. Preprocessing of the data (interpolation) is often necessary to prepare the data for surface-wave mitigation via Fourier methods, and the interpolation itself may be inaccurate due to aliasing.
Another filtering process that is used to reduce the effects of surface waves is adaptive filtering, which attempts to maximize noise suppression, while protecting the desired signal by allowing the algorithm to perform filtering based on an on-the-fly analysis of the data in each record. For example, see Ozbeck, et al., U.S. Pat. No. 6,446,008 B1, which is hereby incorporated herein by reference. However, note that in this filtering process the adaptivity in these methods is to the characteristics of the data in each record, namely the noise level and statistics. While these methods are adaptive in the specific filters they apply to the record, they still require an a priori assumption of the characteristics of body waves and surface waves. Thus, it is still necessary to define the propagation characteristics of the desired output response, e.g. body waves, and the propagation characteristics of the undesired output response, e.g. surface waves, in the f-k space. The adaptive methods so described cannot change the definition of body wave and surface wave propagation characteristics because there is no analysis of those propagation characteristics except perhaps in an ad hoc manner of analyzing those characteristics at a few selected locations, followed by interpolation between those locations.
A further filtering process that is used to reduce the effects of surface waves is phase-matched filtering, which is a method of removing the dispersion characteristics of the surface waves by flattening the surface waves in a seismic record. Phase matching also compresses the long and ringy surface-wave waveform in the time domain by removing the frequency-dependent velocity structure of the surface wave. This produces a surface wave that is not only flat but compact in the time-space domain of the seismic record. This compression of the surface wave is very advantageous because it allows small windows to be used over the limited frequency range of the surface wave to remove the surface wave. In an improvement to narrow time windows, Kim, U.S. Pat. No. 5,781,503, which is hereby incorporated herein by reference, teaches the use of a spatial low-pass filter on the time-aligned and compressed surface-wave data.
Despite the value of aligning and compressing the surface waves, and the subsequent spatial low-pass filtering, it is still necessary with phase-matched filtering of any kind to perform an analysis of the dispersion characteristics of the surface waves. These dispersion characteristics, or frequency dependent velocities, are again traditionally analyzed on some representative records from around the survey area, and ad hoc methods are used to interpolate the dispersion characteristics to locations not analyzed.
FIG. 1 depicts a typical process 100 to mitigate surface waves. Various filtering techniques referred to above may be used in the process 100. The process starts with one or more seismic records 101 for a particular region of interest. In block 102, the records 101 are analyzed at very few, selected, locations in the seismic survey. The analysis involves determining the velocities and dispersion curves at the very few, selected locations. The data resulting from block 102 are sparsely sampled surface-wave properties 103. With this data 103, the process then designs filtering criteria to separate surface waves from body waves in block 104. The resulting sparse filtering criteria 105 are then interpolated by the process in block 106 for every location in the record and for every record in the input data 101. The interpolated criteria 107 are used in block 108 by the process for the mitigation of surface waves in the input data 101 to produce data 109 with mitigated surface waves. Note that in other processes, the surface-wave properties are interpolated for every record instead of the filtering criteria being interpolated, but they result in the same inexact knowledge of the surface wave properties and/or the filtering criteria to separate surface waves from body waves.
In earthquake seismology, surface wave signals are cross-correlated between two seismic stations, and seismic tomography is employed to estimate the group velocities of the surface waves between the stations, for example see Shapiro et al., “High-Resolution Surface-Wave Tomography from Ambient Seismic Noise,” Science Vol. 307, pp. 1615-1618, 2005, which is hereby incorporated herein by reference. This approach is time consuming and computationally expensive, since (1) tomography needs to be employed to overcome sparse spatial sampling, and (2) typically N×N cross-correlations are computed to make full use of available sparse data, where N is the number of seismic stations. Furthermore, in earthquake seismology, surface-wave signals are often directly cross-correlated to obtain the group velocities of the surface waves, without converting the surface-waves into analytic signals. Note that this is possible only when the distance between the seismic stations is much larger than the wavelength of the surface waves, and thus is not applicable to exploration seismology.
In US Patent Application 2005/0152220 A1, which is hereby incorporated herein by reference, Kritski and Amundsen disclose the use of cross-correlation of surface waves recorded on nearby traces in a vertical seismic profile (VSP) to obtain the shear-wave velocity in the region between their depth levels. However, this methodology involves attempting to estimate the entire dispersion curve (shear-wave velocity vs. frequency) by first transforming the data into the wavelet transform domain before correlation.
Similarly, U.S. Pat. No. 6,612,398 B1 to Tokimatsu et al. and U.S. Pat. No. 5,035,144 to Aussel, both of which are hereby incorporated herein by reference, teach the use of cross-correlation to obtain the shear-wave velocity vs. frequency for the fundamental mode of a surface wave. These disclosures obtain the full dispersion curve, and hence first perform a frequency transform before correlating. Aussel constructs the envelope function using the analytic signal many times, namely for each individual frequency component. Note that Tokimatsu operates with small areas at ground level, and Aussel is operative for small quantities of data obtained in nondestructive testing experiments with 2 cm thick material samples.
In U.S. Pat. No. 6,266,620 B1, which is hereby incorporated herein by reference, Baeten and Lemenager teach a method of automatically detecting the cone of the surface waves or ground roll waves. It is the object of that method to identify where on the seismic record surface-wave energy resides, in order that an adaptive filter can be applied to the filtering of that energy. Accordingly, that method relies primarily on amplitudes of surface waves being larger than those of other types of signals, and employs certain types of amplitude thresholding, blocking, and merging to identify and isolate the surface-wave noise. While they perform cross-correlation between the adjacent traces to calculate the velocities of the surface waves in the course of their methodology, the cross-correlation is made between the seismic traces without turning the traces into complex analytic signals. This results in velocity estimates that are in between the phase and the group velocities of the surface waves due to their dispersive nature. Also note that while the methods in U.S. Pat. No. 6,266,620 B1 also turns the traces into analytic signals, the magnitude of the complex analytic signal, or the envelope signal, is only used for instantaneous amplitude estimation. This is because the method relies heavily on certain amplitude criteria. Since the goal of the method is to identify the surface wave noise cone, the output is a range of velocities per trace that supposedly confine all the surface wave noise.
In U.S. Pat. No. 5,241,514 A, which is hereby incorporated herein by reference, Ehlers teaches a method of identifying, characterizing, and mapping seismic noise, which includes surface-wave noise. This method rapidly characterizes surface-wave scatterer location and strength, assuming the surface-wave velocity is neither space nor frequency dependent. As part of the beam steering process, Ehlers teaches estimating the strength and direction of scatterers by varying the velocity of surface waves at each azimuth, in order to identify which combination of velocity and location provides the most energetic display of surface-wave energy. However, no allowance is made for the variability of surface wave along a propagation path for a fixed azimuth. As in the radar field, which that method is based on, velocity along a path is fixed and assumed known for each scan. In the presence of velocity heterogeneity, that method has an ambiguity between scatterer location and velocity. Furthermore, the method does not vary the velocity of the surface waves as a function of frequency, and so is not readily applicable to surface waves that exhibit strong dispersion.